Step by step operation for multi node -- an example with two nodes¶

This document completes the single node operation optimisation with a multinode case (2 nodes)

  • If you need to learn the one node model first with associated equation you have to go there : one node model README it gives planing tools and consumption models
  • If you want the more general with on the two node model, including planing, you can go there : two node models README for step by step learning with two nodes (France and Germany).
  • If you want a more realistic case, you can go there : complete 7 node model README a more complete European model to do prospective analysis
  • 1. optimisation problem
  • 2. Case without storage
  • 3. Case with storage

1. optimisation problem ¶

TODO add congestion constraint¶

\begin{align} &\text{Cost function }& &\min_{x} \sum_z \sum_t \sum_i \pi_{iz} x_{itz}\;\;\; & & \pi_{iz} \text{ marginal cost}\\ &\text{Power limit } & &\text{ s.t.} \;\; 0 \leq x_{itz}\leq a_{itz} \bar{x_{iz}} & &\bar{x_{iz}} \text{ installed power, } a_{itz} \text{ availability}\\ &\text{Meet demand } & & \sum_i x_{itz} \geq C_{tz} && C_{tz} \text{ Consumption}\\ &\text{Stock limit } & &\sum_t x_{it}\leq E_i && E_i=\bar{x_i}*N_i \text{ Energy capacity limit}\\ &\text{ramp limit } & &rc^-_i *x_{it}\leq x_{it}-x_{i(t+1)}\leq rc^+_i *x_{it} && rc^+_i rc^-_i\text{ ramp limit}\\ \end{align}
In [1]:
#region importation of modules
import os
import numpy as np
import pandas as pd
import csv

import datetime
import copy
import plotly.graph_objects as go
import matplotlib.pyplot as plt
from sklearn import linear_model
import sys

import os
import sys
if os.path.basename(os.getcwd())=='Operation_optimisation':
    sys.path.append('../../../')


from Models.Basic_France_models.Operation_optimisation.f_operationModels import *
from functions.f_tools import *
from functions.f_graphicalTools import *
from functions.f_consumptionModels import *

## locally defined models
from Models.Basic_France_Germany_models.Operation_optimisation.f_operationModels import *

#endregion

#region Solver and data location definition

if sys.platform != 'win32':
    myhost = os.uname()[1]
else : myhost = ""
if (myhost=="jupyter-sop"):
    ## for https://jupyter-sop.mines-paristech.fr/ users, you need to
    #  (1) run the following to loanch the license server
    if (os.system("/opt/mosek/9.2/tools/platform/linux64x86/bin/lmgrd -c /opt/mosek/9.2/tools/platform/linux64x86/bin/mosek.lic -l lmgrd.log")==0):
        os.system("/opt/mosek/9.2/tools/platform/linux64x86/bin/lmutil lmstat -c 27007@127.0.0.1 -a")
    #  (2) definition of license
    os.environ["MOSEKLM_LICENSE_FILE"] = '@jupyter-sop'

BaseSolverPath='/Users/robin.girard/Documents/Code/Packages/solvers/ampl_macosx64' ### change this to the folder with knitro ampl ...
## in order to obtain more solver see see https://ampl.com/products/solvers/open-source/
## for eduction this site provides also several professional solvers, that are more efficient than e.g. cbc
sys.path.append(BaseSolverPath)
solvers= ['gurobi','knitro','cbc'] # try 'glpk', 'cplex'
solverpath= {}
for solver in solvers : solverpath[solver]=BaseSolverPath+'/'+solver
solver= 'mosek' ## no need for solverpath with mosek.
#endregion
InputFolder='Data/input/'
GraphicalResultsFolder="GraphicalResults/"

InputConsumptionFolder='../Consumption/Data/'
InputProductionFolder='../Production/Data/'
InputOperationFolder='Data/'
InputEcoAndTech = '../Economic_And_Tech_Assumptions/'

1. Case without storage ¶

In [2]:
#region I - Ramp Ctrs multiple area : loading parameters
Zones="FR_DE_GB_ES"
year=2016
Selected_AREAS=["FR","DE"]
Selected_TECHNOLOGIES=['OldNuke', 'CCG','WindOnShore',"curtailment"] #you'll add 'Solar' after #'NewNuke', 'HydroRiver', 'HydroReservoir','WindOnShore', 'WindOffShore', 'Solar', 'Curtailement'}

#### reading CSV files
TechParameters = pd.read_csv(InputOperationFolder+'Gestion_MultiNode_DE-FR_AREAS_TECHNOLOGIES.csv',sep=',',decimal='.',comment="#",skiprows=0).set_index(["AREAS","TECHNOLOGIES"])
areaConsumption = pd.read_csv(InputConsumptionFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
                                sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["AREAS","Date"])
availabilityFactor = pd.read_csv(InputProductionFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
                                sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["AREAS","Date","TECHNOLOGIES"])

ExchangeParameters = pd.read_csv(InputEcoAndTech+'Hypothese_DE-FR_AREAS_AREAS.csv',sep=',',decimal='.',skiprows=0,comment="#").set_index(["AREAS","AREAS.1"])
#ExchangeParameters.loc[ExchangeParameters.AREAS=="FR",'maxExchangeCapacity']=90000 ## margin to make everything work
#ExchangeParameters.loc[ExchangeParameters.AREAS=="DE",'maxExchangeCapacity']=90000 ## margin to make everything work
#### Selection of subset
TechParameters=TechParameters.loc[(Selected_AREAS,Selected_TECHNOLOGIES),:]
areaConsumption=areaConsumption.loc[(Selected_AREAS,slice(None)),:]
availabilityFactor=availabilityFactor.loc[(Selected_AREAS,slice(None),Selected_TECHNOLOGIES),:]
TechParameters.loc[(slice(None),'CCG'),'capacity']=100000 ## margin to make everything work
TechParameters.loc[(slice(None),"OldNuke"),'RampConstraintMoins']=0.001 ## a bit strong to put in light the effect
TechParameters.loc[(slice(None),"OldNuke"),'RampConstraintPlus']=0.001 ## a bit strong to put in light the effect
#endregion

/var/folders/zt/7cz1y9md79l_h08bbqymt4w9z8xlw7/T/ipykernel_42278/1973306466.py:20: FutureWarning: The behavior of indexing on a MultiIndex with a nested sequence of labels is deprecated and will change in a future version. `series.loc[label, sequence]` will raise if any members of 'sequence' or not present in the index's second level. To retain the old behavior, use `series.index.isin(sequence, level=1)`
  availabilityFactor=availabilityFactor.loc[(Selected_AREAS,slice(None),Selected_TECHNOLOGIES),:]
In [3]:
#region I - Ramp Ctrs multiple area : solving and loading results
### small data cleaning
availabilityFactor.availabilityFactor[availabilityFactor.availabilityFactor>1]=1
model = GetElectricSystemModel_GestionMultiNode(Parameters = {"areaConsumption" : areaConsumption,
                                                               "availabilityFactor" : availabilityFactor,
                                                               "TechParameters" : TechParameters,
                                                              "ExchangeParameters":ExchangeParameters})
opt = SolverFactory(solver)
results=opt.solve(model)
Variables=getVariables_panda_indexed(model)
production_df=EnergyAndExchange2Prod(Variables)

### Check sum Prod = Consumption
Delta= production_df.sum(axis=1)-areaConsumption.areaConsumption
abs(Delta).sum()

print(production_df.loc[('DE',slice(None)),'OldNuke'].diff(1).max()/TechParameters.loc[("DE","OldNuke"),"capacity"])
print(production_df.loc[('DE',slice(None)),'OldNuke'].diff(1).min()/TechParameters.loc[("DE","OldNuke"),"capacity"])
print(production_df.loc[('FR',slice(None)),'OldNuke'].diff(1).max()/TechParameters.loc[("FR","OldNuke"),"capacity"])
print(production_df.loc[('FR',slice(None)),'OldNuke'].diff(1).min()/TechParameters.loc[("FR","OldNuke"),"capacity"])
## adding dates

fig=MyAreaStackedPlot(production_df,Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename=GraphicalResultsFolder+'file.html') ## offline
fig.show()

production_df.sum(axis=0)/10**6 ### energies produites TWh
production_df.groupby(by="AREAS").sum()/10**6 ### energies produites TWh
production_df.groupby(by=["AREAS"]).max()

Constraints= getConstraintsDual_panda(model)
Constraints.keys()
Constraints['energyCtr']
#endregion

0.0009468354430380174
-0.0009495253164557114
0.0009468354430380656
-0.0009495253164557017
Out[3]:
AREAS Date energyCtr
0 DE 2016-06-29 20:00:00 68.8
1 DE 2016-09-02 01:00:00 68.8
2 DE 2016-12-22 17:00:00 68.8
3 DE 2016-06-26 04:00:00 68.8
4 DE 2016-09-07 19:00:00 68.8
... ... ... ...
17563 FR 2016-07-11 06:00:00 68.8
17564 FR 2016-01-28 06:00:00 0.0
17565 FR 2016-06-05 08:00:00 68.8
17566 FR 2016-07-16 23:00:00 68.8
17567 FR 2016-04-02 10:00:00 68.8

17568 rows × 3 columns

2. Case with storage ¶

In [4]:
#region II Ramp+Storage Multi area : loading parameters
Zones="FR_DE_GB_ES"
year=2016
Selected_AREAS=["FR","DE"]
Selected_TECHNOLOGIES=['OldNuke', 'CCG','WindOnShore',"curtailment"] #you'll add 'Solar' after #'NewNuke', 'HydroRiver', 'HydroReservoir','WindOnShore', 'WindOffShore', 'Solar', 'Curtailement'}

#### reading CSV files
TechParameters = pd.read_csv(InputOperationFolder+'Gestion_MultiNode_DE-FR_AREAS_TECHNOLOGIES.csv',
                             sep=',',decimal='.',comment="#",skiprows=0).set_index(["AREAS","TECHNOLOGIES"])
StorageParameters = pd.read_csv(InputOperationFolder+'Gestion_MultiNode_AREAS_DE-FR_STOCK_TECHNO.csv',sep=',',decimal='.',comment="#",skiprows=0).set_index(["AREAS","STOCK_TECHNO"])
areaConsumption = pd.read_csv(InputConsumptionFolder+'areaConsumption'+str(year)+'_'+str(Zones)+'.csv',
                                sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["AREAS","Date"])
availabilityFactor = pd.read_csv(InputProductionFolder+'availabilityFactor'+str(year)+'_'+str(Zones)+'.csv',
                                sep=',',decimal='.',skiprows=0,parse_dates=['Date']).set_index(["AREAS","Date","TECHNOLOGIES"])

ExchangeParameters = pd.read_csv(InputEcoAndTech+'Hypothese_DE-FR_AREAS_AREAS.csv',sep=',',decimal='.',skiprows=0,comment="#").set_index(["AREAS","AREAS.1"])
#ExchangeParameters.loc[ExchangeParameters.AREAS=="FR",'maxExchangeCapacity']=90000 ## margin to make everything work
#ExchangeParameters.loc[ExchangeParameters.AREAS=="DE",'maxExchangeCapacity']=90000 ## margin to make everything work
#### Selection of subset
TechParameters=TechParameters.loc[(Selected_AREAS,Selected_TECHNOLOGIES),:]
areaConsumption=areaConsumption.loc[(Selected_AREAS,slice(None)),:]
availabilityFactor=availabilityFactor.loc[(Selected_AREAS,slice(None),Selected_TECHNOLOGIES),:]
TechParameters.loc[(slice(None),'CCG'),'capacity']=100000 ## margin to make everything work
TechParameters.loc[(slice(None),"OldNuke"),'RampConstraintMoins']=0.01 ## a bit strong to put in light the effect
TechParameters.loc[(slice(None),"OldNuke"),'RampConstraintPlus']=0.02 ## a bit strong to put in light the effect
#endregion

/var/folders/zt/7cz1y9md79l_h08bbqymt4w9z8xlw7/T/ipykernel_42278/3663365949.py:22: FutureWarning:

The behavior of indexing on a MultiIndex with a nested sequence of labels is deprecated and will change in a future version. `series.loc[label, sequence]` will raise if any members of 'sequence' or not present in the index's second level. To retain the old behavior, use `series.index.isin(sequence, level=1)`
In [5]:
#region II Ramp+Storage multi area : solving and loading results
model= GetElectricSystemModel_GestionMultiNode_withStorage(Parameters = {"areaConsumption" : areaConsumption,
                                                               "availabilityFactor" : availabilityFactor,
                                                               "TechParameters" : TechParameters,
                                                              "ExchangeParameters":ExchangeParameters,
                                                              "StorageParameters": StorageParameters})
if solver in solverpath :  opt = SolverFactory(solver,executable=solverpath[solver])
else : opt = SolverFactory(solver)
results=opt.solve(model)
Variables = getVariables_panda_indexed(model)
Constraints = getConstraintsDual_panda(model)

production_df=EnergyAndExchange2Prod(Variables)
stockage=Variables['storageOut'].pivot(index=['AREAS','Date'],columns='STOCK_TECHNO',values='storageOut').sum(axis=1)-Variables['storageIn'].pivot(index=['AREAS','Date'],columns='STOCK_TECHNO',values='storageIn').sum(axis=1)
areaConsumption['Storage']=stockage
areaConsumption['NewConsumption']=areaConsumption['areaConsumption']+stockage

### Check sum Prod = Consumption
Delta=(production_df.sum(axis=1) - areaConsumption.NewConsumption); ## comparaison à la conso incluant le stockage
abs(Delta).max()
production_df.loc[:,'Storage'] = -areaConsumption["Storage"] #### ajout du stockage comme production
Delta=(production_df.sum(axis=1) - areaConsumption.areaConsumption);
abs(Delta).max()

fig=MyAreaStackedPlot(production_df,Conso=areaConsumption)
fig=fig.update_layout(title_text="Production électrique (en KWh)", xaxis_title="heures de l'année")
#plotly.offline.plot(fig, filename=GraphicalResultsFolder+'file.html') ## offline
fig.show()

abs(areaConsumption["Storage"]).groupby(by="AREAS").sum() ## stockage
production_df.groupby(by="AREAS").sum()/10**6 ### energies produites TWh
production_df[production_df>0].groupby(by="AREAS").sum()/10**6 ### energies produites TWh
production_df.groupby(by="AREAS").max()/1000 ### Pmax en GW ### le stockage ne fait rien en Allemagne ??? bizarre
production_df.groupby(by="AREAS").min()/1000 ### Pmax en GW
#endregion
Out[5]:
CCG OldNuke WindOnShore curtailment DE FR Storage
AREAS
DE 0.000000e+00 4.310174 0.167386 0.0 0.0 -9.0 -5.0
FR -9.094947e-16 33.523576 0.070911 0.0 -9.0 0.0 -5.0